
rm(list = ls())

# Load required libraries
library(haven)
library(dplyr)
library(ggplot2)
library(lmtest)
library(sandwich)
library(stargazer)
library(tidyverse)
library(broom)
library(RDHonest)

# Function to cluster standard errors
cluster_se <- function(model, data, type = "HC1") {
  idx <- as.integer(rownames(model.frame(model)))
  vc  <- sandwich::vcovCL(model, cluster = data$cluster[idx], type = type)
  lmtest::coeftest(model, vcov. = vc)
}

# Function to run placebo RD models and output LaTeX table
run_placebo_models <- function(data_path, treat_var, dosage_var, label, col_labels) {
  cep_all <- read_dta(data_path)
  dosage_ranges <- 1:5
  
  # Create list of samples for each window
  sample_list <- lapply(dosage_ranges, function(i) {
    subset(cep_all, cep_all[[dosage_var]] >= -i & cep_all[[dosage_var]] <= i)
  })
  
  # Fit models
  models_lm <- lapply(seq_along(sample_list), function(i) {
    if (i == 1) {
      lm(as.formula(paste0("left ~ ", treat_var, " + as.factor(encuesta)")), data = sample_list[[i]])
    } else {
      lm(as.formula(paste0("left ~ ", treat_var, " + ", dosage_var, " + ",
                           dosage_var, ":", treat_var, " + as.factor(encuesta)")),
         data = sample_list[[i]])
    }
  })
  
  # Clustered SEs
  models_cl <- Map(cluster_se, models_lm, sample_list)
  
  # Summary stats
  sample_sizes <- sapply(models_lm, nobs)
  r_squared <- sapply(models_lm, function(m) summary(m)$r.squared)
  mean_dv <- sapply(models_lm, function(m) mean(m$model$left, na.rm = TRUE))
  
  # Stargazer output
  stargazer(models_cl,
            type = "latex",
            title = paste("Table G1: RD Estimates 1973 Coup on Left-wing Ideological Identification (Placebo cutoff", label, ")"),
            keep = c("^treat_placebo_1$", "^treat_placebo1$"),
            column.labels = col_labels,
            covariate.labels = "Coup",
            dep.var.caption = "Dependent Variable: Left Identification",
            dep.var.labels.include = FALSE,
            add.lines = list(
              c("Obs.", sample_sizes),
              c("$R^2$", round(r_squared, 3)),
              c("Mean DV", round(mean_dv, 3))
            ),
            align = TRUE)
}

# Run for Placebo Cutoff +1
run_placebo_models(
  data_path = "Data/cep_all_placebo_cut.dta",
  treat_var = "treat_placebo1",
  dosage_var = "treatment_dosage_coup_placebo1",
  label = "+1",
  col_labels = c("55-57", "54-58", "59-53", "60-52", "51-61")
)

# Run for Placebo Cutoff -1
run_placebo_models(
  data_path = "Data/cep_all_placebo_cut.dta",
  treat_var = "treat_placebo_1",
  dosage_var = "treatment_dosage_coup_placebo_1",
  label = "-1",
  col_labels = c("53-55", "52-56", "51-57", "50-58", "49-59")
)

print('Table G1 complete')

